c     newestad.f  :  obtiene informacion de los snapshots de Gadget2
c
      implicit none
      integer nmax
      parameter (nmax=1000000)
      integer isnap,igrid,id(nmax),npart(6),npar(6),ntot,sf,fb
      integer cool,nfiles,i,j,k,imaxj
      real*8 twopi, k2, mjup, pi, masaprin, masa2
      real*8 time,z,box,om,la,h,mass(6)
      real*8 axj,exj,incj,zmj,omegaj,omj, posj(6)
      real pos(1:3,nmax),vel(1:3,nmax),masadw(nmax),u(nmax)
      real sl(nmax),densi(nmax),disco(6),masdisk,comp(3)
      real estad(1:11,150),nestad(1:6,150),elemt(3,nmax)
      real pos1(1:3,nmax),vel1(1:3,nmax),r1(nmax),v1(nmax)
      real ai, af, da, ei, ef, de, wi, wf, dw, mastot,li,lf,dl
      character(len=50) arch_big,arch_in,arch_out01,arch_out02
      character(len=3)  next
      character(len=4)  next1
      character(len=5)  next2
      byte c5(96)
c
      common /gen/ twopi, k2, mjup, pi
c
c      apertura del archivo de evolucion de la compañera, del disco y de la perdida de masa
      open (44,file='evolu.dat',status='replace')
      close (44)
c
c        inicio de la lectura de los snapshot, apertura de archivos
        isnap=60000
       do igrid=0,isnap  ! gran do
       open (55,status='scratch')
       if (igrid.lt.10) write (55,'(3i1)') 0,0,igrid
       if (igrid.ge.10.and.igrid.lt.100) write (55,'(i1,i2)') 0,igrid
       if (igrid.ge.100.and.igrid.lt.1000) write (55,'(i3)') igrid
       if (igrid.ge.1000.and.igrid.lt.10000) write (55,'(i4)') igrid
       if (igrid.ge.10000) write (55,'(i5)') igrid
       rewind (55)
       if (igrid.ge.1000) then
        if (igrid.ge.10000) then
        read (55,'(a5)') next2
        else
        read (55,'(a4)') next1
        end if
       else
       read (55,'(a3)') next
       end if
       close (55)
       if (igrid.ge.1000) then
        if (igrid.ge.10000) then
        arch_big = 'snapshot2_'//next2//''
        else
        arch_big = 'snapshot2_'//next1//''
        end if
       else
       arch_big = 'snapshot2_'//next//''
       end if
      arch_in = arch_big
c
c      apertura de los archivos con informacion. histo --> posee los histogramas de a, e y w de las particulas del disco (con a<10)
c                                                general --> da distribuciones de masa, densidad, a, e, w y Q
       if (igrid.ge.1000) then
        if (igrid.ge.10000) then
      arch_out01 = 'histo_'//next2//'.dat'
      arch_out02 = 'general_'//next2//'.dat'
        else
      arch_out01 = 'histo_'//next1//'.dat'
      arch_out02 = 'general_'//next1//'.dat'
        end if
       else
      arch_out01 = 'histo_'//next//'.dat'
      arch_out02 = 'general_'//next//'.dat'
       end if
c
      open (1,file=arch_in,form='unformatted',status='old')
      open (31,file=arch_out01,status='replace')
      open (32,file=arch_out02,status='replace')
c
c     constantes y parametros globales.
      twopi = 8.0*datan(1.0d0)
      pi = twopi/2.0d0
      k2    = 0.01702029895**2
      mjup  = 9.45d-4
c
c     lee header de cada snapshot.
      read (1) npar,mass,time,z,sf,fb,npart,cool,nfiles,box,om,la,h,c5
c
c     calcula el numero total de cuerpos.
      ntot = sum(npar)
c
c     lee datos de las particulas.
      read (1) pos(1:3,1:ntot)
      read (1) vel(1:3,1:ntot)
      read (1) id(1:ntot)
      read (1) masadw(1:ntot)
      read (1) u(1:npar(1))
      if (igrid.lt.1) mastot=masadw(20)*npar(1) ! masa total del disco
c     deja en 0 algunas variables a utilizar
      disco=0.0
      estad=0.0
      nestad=0.0
      masdisk=0.0
c     busca la estrella principal
      do i=1,ntot
      if (masadw(i).gt.1.5) imaxj=i
      end do
c      Inicia el analisis de los datos
      do i=1,ntot
c
c     Paso a sistema astorcentrico de m1 y m2 y calculo de modulos de posicion y velocidad
      call cambio(imaxj,pos,vel,pos1,vel1,ntot,nmax,r1,v1,i) ! posicion y velocidad original,posicion y velocidad astrocentrica 
c      call modulos(pos1,vel1,ntot,nmax,r1,v1)  ! calcula el modulo del vector posicion y velocidad de m1
c
c       determino los elementos orbitales de las particulas, carga los histogramas y los elementos orbitales de la compañera.
        if (masadw(i).lt.1.5) then
         posj(1)=pos1(1,i)
         posj(2)=pos1(2,i)
         posj(3)=pos1(3,i)
         posj(4)=vel1(1,i)
         posj(5)=vel1(2,i)
         posj(6)=vel1(3,i)
         masaprin=1.6d0
         masa2=masadw(i)
         call elem(masaprin,masa2,posj,axj,exj,incj,
     *    zmj,omegaj,omj) ! calcula los elementos orbitales de cada particula
         elemt(1,i)=axj
         elemt(2,i)=exj
         elemt(3,i)=omegaj
c
          if (masadw(i).gt.0.2)then
          comp(1)=axj
          comp(2)=exj
          comp(3)=omegaj
          end if
c
          if (masadw(i).lt.0.2)then
           if (r1(i).lt.10.0) then
           masdisk=masdisk+masadw(i)
           disco(1)=axj+disco(1)
           disco(2)=exj+disco(2)
           disco(3)=omegaj+disco(3)
           end if
          end if
        end if
c
c       calcula el perfil de masa, densidad y la media por intervalo de r       
        if (masadw(i).lt.0.1) then
         li=0.0001
         do j=1,101
         lf=li+0.1
         dl=0.1
         estad(1,j)=(lf+li)/2.0
         if (r1(i).ge.li.and.r1(i).lt.lf) then
         estad(2,j)=estad(2,j)+1.0
         estad(3,j)=estad(3,j)+masadw(i)
         estad(4,j)=estad(4,j)+elemt(1,i)
         estad(5,j)=estad(5,j)+elemt(2,i)
         estad(6,j)=estad(6,j)+elemt(3,i)
         estad(7,j)=estad(7,j)+masadw(i)/twopi/estad(1,j)/dl
         end if
         li=li+0.1
         end do
        end if
c
c       calcula la media de a, e y w del disco completo (para r<10)
        if (masadw(i).lt.0.1) then
         ai=0.0
         ei=0.0
         wi=0.0
         do k=1,101
         af=ai+0.1
         da=0.1
         nestad(1,k)=(af+ai)/2.0
         ef=ei+0.01
         de=0.01
         nestad(3,k)=(ef+ei)/2.0
         wf=wi+0.01*twopi
         dw=0.01*twopi
         nestad(5,k)=(wf+wi)/2.0
         if (elemt(1,i).ge.ai.and.elemt(1,i).lt.af) then
         nestad(2,k)=nestad(2,k)+1.0
         end if
         ai=ai+0.1
         if (elemt(2,i).ge.ei.and.elemt(2,i).lt.ef) then
         nestad(4,k)=nestad(4,k)+1.0
         end if
         ei=ei+0.01
         if (elemt(3,i).ge.wi.and.elemt(3,i).lt.wf) then
         nestad(6,k)=nestad(6,k)+1.0
         end if
         wi=wi+0.01*twopi
         end do
c
c
        end if
c
c
      end do
c
c         completa el calculo de la media para el disco completo
          disco(1)=disco(1)/(ntot-2)
          disco(2)=disco(2)/(ntot-2)
          disco(3)=disco(3)/(ntot-2)
c      
c       escribe el archivo de histogramas
        do k=1,101
        write (31,100) nestad(1,k),nestad(2,k),nestad(3,k),
     *   nestad(4,k),nestad(5,k),nestad(6,k) ! archivo de a vs N, e vs N y w vs N (este archivo es uno por snapshot)
        end do 
c
c      completa el calculo de la media por intervalo de r y calcula Q
        do j=1,101  !  calculo de la media de cada intervalo y  Q
        if (estad(2,j).gt.1.0) then
         estad(4,j)=estad(4,j)/estad(2,j)
         estad(5,j)=estad(5,j)/estad(2,j)
         estad(6,j)=estad(6,j)/estad(2,j)
         estad(8,j)=0.05*1.6/pi/estad(1,j)/estad(1,j)/estad(7,j)
        else
         estad(4,j)=0.0
         estad(5,j)=0.0
         estad(6,j)=0.0
         estad(8,j)=0.0
        end if
        end do 
c
c nuevo do total para calcular desviaciones
      do i=1,ntot
c      desviaciones para el disco completo
        if (masadw(i).lt.0.2)then
          if (r1(i).lt.10.0) then
          disco(4)=disco(4)+(elemt(1,i)-disco(1))**2
          disco(5)=disco(5)+(elemt(2,i)-disco(2))**2
          disco(6)=disco(6)+(elemt(3,i)-disco(3))**2
          end if
        end if
c       calculo de desviaciones para intervalos de r
        if (masadw(i).lt.0.1) then
         li=0.0001
         do j=1,101
         lf=li+0.1
         dl=0.1
         if (r1(i).ge.li.and.r1(i).lt.lf) then
         if (estad(2,j).gt.1.0) then
         estad(9,j)=estad(9,j)+(elemt(1,i)-estad(4,j))**2
         estad(10,j)=estad(10,j)+(elemt(2,i)-estad(5,j))**2
         estad(11,j)=estad(11,j)+(elemt(3,i)-estad(6,j))**2
         else
         estad(9,j)=0.0
         estad(10,j)=0.0
         estad(11,j)=0.0
         end if
         end if
         li=li+0.1
         end do
        end if
      end do
c
c         completa el calculo de desviaciones para el disco total
          disco(4)=sqrt(disco(4)/(ntot-3))  ! desviacion estandar para el disco total 
          disco(5)=sqrt(disco(5)/(ntot-3))
          disco(6)=sqrt(disco(6)/(ntot-3))
c       escritura del archivo de evolucion, y escritura en pantalla para llevar control.
      open (44,file='evolu.dat',status='old',access='append')
            write (44,100) time,comp(1),comp(2),comp(3),
     *              disco(1),disco(2),disco(3),disco(4),
     *        disco(5),disco(6),masdisk/mastot,time/56.25 !archivo de evolucion
            write (*,100) time,comp(1),comp(2),comp(3),
     *              disco(1),disco(2),disco(3),disco(4),
     *              disco(5),disco(6),masdisk/mastot
             close (44)
c
c      completa el calculo de desviaciones par intervalos de r
        do j=1,101  !  calculo de la desv standar de cada intervalo
        if (estad(2,j).gt.1.0) then
        estad(9,j)=sqrt(estad(9,j)/(estad(2,j)-1.0))
        estad(10,j)=sqrt(estad(10,j)/(estad(2,j)-1.0))
        estad(11,j)=sqrt(estad(11,j)/(estad(2,j)-1.0))
        else
        estad(9,j)=0.0
        estad(10,j)=0.0
        estad(11,j)=0.0
        end if
        end do 
c
c        escritura en el archivo general de los perfiles de m, sigma, a, e, w y Q
        do j=1,101  ! escritura del archivo ppal
       write (32,100) estad(1,j),estad(2,j),estad(3,j),
     *     estad(4,j),estad(5,j),estad(6,j),estad(7,j),
     *    estad(8,j),estad(9,j),estad(10,j),estad(11,j) ! punto medio, N, masa, a,e, w, sigma, Q, dsa, dse, dsw
        end do
c
c
       end do    ! final del lazo de snapshot
c
c      formato
 100  format (1p20d15.5)
c
      end


      subroutine cambio(imax,zpos,zvel,zpos1,zvel1,ntotz,nmaxz,
     *   rtz,vtz,ij)
      integer ntotz,ij,imax,nmaxz
      real zpos(1:3,ntotz),zvel(1:3,nmaxz)
      real rtz(nmaxz),vtz(nmaxz)
      real zpos1(1:3,nmaxz),zvel1(1:3,nmaxz)
c
      common /gen/ twopi, k2, mjup, pi
c
c       do ij = 1,ntotz
       zpos1(1,ij)=zpos(1,ij)-zpos(1,imax)
       zpos1(2,ij)=zpos(2,ij)-zpos(2,imax)
       zpos1(3,ij)=zpos(3,ij)-zpos(3,imax)
       zvel1(1,ij)=zvel(1,ij)-zvel(1,imax)
       zvel1(2,ij)=zvel(2,ij)-zvel(2,imax)
       zvel1(3,ij)=zvel(3,ij)-zvel(3,imax)
              rtz(ij)=sqrt(zpos1(1,ij)**2+
     *      zpos1(2,ij)**2+zpos1(3,ij)**2)
       vtz(ij)=sqrt(zvel1(1,ij)**2+
     *      zvel1(2,ij)**2+zvel1(3,ij)**2)
c       end do
c       end do
c
       end


      subroutine modulos(zpos1,zvel1,ntotz,nmaxz,rtz,vtz)
      integer ij,ntotz,nmaxz
      real rtz(nmaxz),vtz(nmaxz)
      real zpos1(1:3,nmaxz),zvel1(1:3,nmaxz)
c
      common /gen/ twopi, k2, mjup, pi
c
       do ij = 1,ntotz
       rtz(ij)=sqrt(zpos1(1,ij)**2+
     *      zpos1(2,ij)**2+zpos1(3,ij)**2)
       vtz(ij)=sqrt(zvel1(1,ij)**2+
     *      zvel1(2,ij)**2+zvel1(3,ij)**2)
       end do
c
       end


      subroutine elem (m1,m2,xc,a,e,inc,capm,omega,capom)
c
      implicit real*8 (a-h,k-z)
      parameter (pi = 3.14159265358979d0)
      parameter (piby2 = pi/2.0d0)
      parameter (degrad = 180.0d0 / pi)
c      real xc(6)
      real*8 xc(6),inc
c      common /gen/ twopi,cero,uno,dos,tres2,uno2,g,rad,error,unoc2,gsum0
c
      error = 1.0d-13 
      twopi = 8.0d0*datan(1.0d0)
      cero  = 0.0d0
      uno   = 1.0d0
      dos   = 2.0d0
      uno2  = uno/dos
      tres2 = 3.0d0/2.0d0
c      pi    = uno2*twopi
      g     = dsqrt(39.477d0) !1.720209895d-02
      rad   = twopi/360.0
c
      msum=m1+m2
      gmsum = g*g*msum
      x  = xc(1)
      y  = xc(2)
      z  = xc(3)
      vx = xc(4)
      vy = xc(5)
      vz = xc(6)
      
c     compute the angular momentum h, and thereby the inclination inc.
      hx = y*vz - z*vy
      hy = z*vx - x*vz
      hz = x*vy - y*vx
      h2 = hx*hx + hy*hy +hz*hz
      h  = sqrt(h2)
      inc = acos(hz/h)
      
c     compute longitude of ascending node capom and the argument of latitude u.
      fac = sqrt(hx**2 + hy**2)/h
      if (fac.lt.error) then
       capom = 0.d0
       u = atan2(y,x)
       if(abs(inc - pi).lt.10.d0*error) u = -u
      else
       capom = atan2(hx,-hy)	  
       u = atan2 ( z/sin(inc) , x*cos(capom) + y*sin(capom))
      endif
      if(capom .lt. 0.d0) capom = capom + 2.d0*pi
      if(u .lt. 0.d0) u = u + 2.d0*pi
      
c     compute the radius r, vel. squared v2, the dot product rdotv & energy.
      r = sqrt(x*x + y*y + z*z)
      v2 = vx*vx + vy*vy + vz*vz
      v = sqrt(v2)
      vdotr = x*vx + y*vy + z*vz
      energy = 0.5d0*v2 - gmsum/r
      
c     determine type of conic section and label it via ialpha
      if(abs(energy*r/gmsum) .lt. sqrt(error)) then
       ialpha = 0
      else
       if(energy .lt. 0.d0) ialpha = -1 
       if(energy .gt. 0.d0) ialpha = +1
      endif

c     depending on the conic type, determine the remaining elements
c     ellipse :
      if(ialpha .eq. -1) then
       a = -0.5d0*gmsum/energy  
       fac = 1.d0 - h2/(gmsum*a)
       if (fac .gt. error) then
        e = sqrt ( fac )
        face =(a-r)/(a*e)
c     apr. 16/93 : watch for case where face is slightly outside unity
        if ( face .gt. 1.d0) then
         cape = 0.d0
        else
         if ( face .gt. -1.d0) then
          cape = acos( face )
         else
          cape = pi
         endif
        endif
        if ( vdotr .lt. 0.d0 ) cape = 2.d0*pi - cape
        cw = (cos( cape) -e)/(1.d0 - e*cos(cape))
        sw = sqrt(1.d0 - e*e)*sin(cape)/(1.d0 - e*cos(cape))
        w = atan2(sw,cw)
        if(w .lt. 0.d0) w = w + 2.d0*pi
       else
        e = 0.d0
        w = u
        cape = u
       endif
       capm = cape - e*sin (cape)
       omega = u - w
       if(omega .lt. 0.d0) omega = omega + 2.d0*pi
       omega = omega - int(omega/(2.d0*pi))*2.d0*pi 	 
      endif
   
c     hyperbola
      if(ialpha .eq. +1) then
       a = +0.5d0*gmsum/energy  
       fac = h2/(gmsum*a)
       if (fac .gt. error) then
        e = sqrt ( 1.d0 + fac )
        tmpf = (a+r)/(a*e)
        if(tmpf.lt.1.0d0) then
         tmpf = 1.0d0
        endif
        capf = log(tmpf + sqrt(tmpf*tmpf -1.d0))
        if ( vdotr .lt. 0.d0 ) capf = - capf
        cw = (e - cosh(capf))/(e*cosh(capf) - 1.d0 )
        sw = sqrt(e*e - 1.d0)*sinh(capf)/(e*cosh(capf) - 1.d0 )
        w = atan2(sw,cw)
        if(w .lt. 0.d0) w = w + 2.d0*pi
       else
c     we only get here if a hyperbola is essentially a parabola
c     so we calculate e and w accordingly to avoid singularities
        e = 1.d0
        tmpf = 0.5d0*h2/gmsum
        w = acos(2.d0*tmpf/r -1.d0)
        if ( vdotr .lt. 0.d0) w = 2.d0*pi - w
        tmpf = (a+r)/(a*e)
        capf = log(tmpf + sqrt(tmpf*tmpf -1.d0))
       endif
       capm = e * sinh(capf) - capf
       omega = u - w
       if(omega .lt. 0.d0) omega = omega + 2.d0*pi
       omega = omega - int(omega/(2.d0*pi))*2.d0*pi 	 
      endif
  
c     parabola : ( note - in this case we use "a" to mean pericentric distance)
      if(ialpha .eq. 0) then
       a =  0.5d0*h2/gmsum  
       e = 1.d0
       w = acos(2.d0*a/r -1.d0)
       if ( vdotr .lt. 0.d0) w = 2.d0*pi - w
       tmpf = tan(0.5d0 * w)
       capm = tmpf* (1.d0 + tmpf*tmpf/3.d0)
       omega = u - w
       if(omega .lt. 0.d0) omega = omega + 2.d0*pi
       omega = omega - int(omega/(2.d0*pi))*2.d0*pi 	 
      endif
c     
      return
      end


